source("/Users/benphillips/evo-dispersal/Spinifex/Space Grass Sim 20100422.R")

##  Initialise the model  ##
ngens<-100
max.pop<-1000000 # total maximum population size for all species combined.  Determines grid-cell carrying capacity
init.or<-c(0,0,0,0,10,0,0)
init.os<-c(1600, 0,0,0,10,0,0)
init.fr<-c(1200, 0,0,0,0,0,10)
lh.or<-matrix(c(0,0,0,70,0,0,0,0,0,0,0,0,0,0,0.58,0,0,0,0,0,0,0,0.9,0.9,0,0,0,0,0,0,0,0.9,0.90,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), nrow=7, byrow=T)
lh.os<-matrix(c(0,200,400,800,1600,0,0,0,0,0,0,0,0,0,0, 0.9,0,0,0,0,0,0,0,0.9,0,0,0,0,0,0,0,0.9,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0), nrow=7, byrow=T)
lh.fr<-matrix(c(0,0,0,150,300,600,1200,0,0,0,0,0,0,0,0,0,0,0,0,0,0, 0,0.9,0.9,0,0,0,0,0,0,0,0.9,0,0,0,0,0,0,0,0.9,0,0,0,0,0,0,0,0.9,0.9), nrow=7, byrow=T)
cl.size<-1 #cell size: must divide into space.size perfectly. One is a sensible default!
sp.size<-10 #Space size
f.scale<-5 #fire scale
cond.prob.fire<-0.05
steps<-sp.size/cl.size # Number of cells along each edge of the "square" space (= periodicity of toroidal space).
K<-floor(max.pop/(steps^2)) # per grid-cell carrying capacity
disp.prob.or<-0.05 #dispersal probability
dp.seeder<-0.005


############################

### Run the model ###
popsize.or<-c()
popsize.os<-c()
popsize.fr<-c()
k.prop.occ<-c()
neigh.grid<-neighbours.init(sp.size, cl.size)
for(k in 1:11){ 
  print(k) 
  kpop.or<-0
  kpop.os<-0
  kpop.fr<-0	
  rep.prop.occ<-c()
  for(j in 1:20){
    rep.pop.or<-c()
    rep.pop.os<-c()
    rep.pop.fr<-c()
    grid<-space.init(sp.size, cl.size, init.or, init.os, init.fr)
    f.grid<-space.fire.init(sp.size, cl.size) # fire pattern for that year
    for(i in 1:ngens){
      f.grid<-fire.grid(f.grid, (k-1)/10, cond.prob.fire, f.scale, cl.size, sp.size)
      for (x in 1:steps){
        for (y in 1:steps){
          # reproduction
          grid[,,,x,y]<-cell.update(grid[,,1,x,y], grid[,,2,x,y], grid[,,3,x,y], lh.or, lh.os, lh.fr, fire.grid.cell=f.grid[x,y])
        }	
      }
      # Dispersal 
      grid<-comb.disp(grid, neigh.grid, disp.prob.or, dp.seeder)
      # collect popsizes in last five generations
      if((ngens-i)<=5) {
        rep.pop.or<-c(rep.pop.or, sum(grid[3:5,1,1,,]))
        rep.pop.os<-c(rep.pop.os, sum(grid[2:5,1,2,,]))
        rep.pop.fr<-c(rep.pop.fr, sum(grid[3:7,1,3,,]))
      }
      if (i==ngens) prop.occ<-prop.occ.collect(grid, steps)
    }
    rep.prop.occ<-rbind(rep.prop.occ, prop.occ)
    rep.pop.or<-mean(rep.pop.or)
    rep.pop.os<-mean(rep.pop.os)
    rep.pop.fr<-mean(rep.pop.fr)	
    kpop.or<-c(kpop.or, rep.pop.or)
    kpop.os<-c(kpop.os, rep.pop.os)
    kpop.fr<-c(kpop.fr, rep.pop.fr)	
  }
  rep.prop.occ<-apply(rep.prop.occ, 2, mean)
  k.prop.occ<-rbind(k.prop.occ, rep.prop.occ)
  popsize.or<-cbind(popsize.or, kpop.or)
  popsize.os<-cbind(popsize.os, kpop.os)
  popsize.fr<-cbind(popsize.fr, kpop.fr)
}
summ.or<-apply(popsize.or, 2, mean)
summ.os<-apply(popsize.os, 2, mean)
summ.fr<-apply(popsize.fr, 2, mean)



par(cex.lab=1.5, mfrow=c(2,1), mar=c(5, 5, 4, 2))
plot(seq(0, 1, 0.1), summ.or, ylim=c(0,max(c(summ.or, summ.os, summ.fr))), xlab="Conditional probability of a fire", ylab="Mean population size", bty="l")
points(seq(0, 1, 0.1), summ.os, col="red")
points(seq(0, 1, 0.1), summ.fr, col="blue")
legend(0.5, max(c(summ.or, summ.os, summ.fr)), legend=c("obligate resprouter", "obligate seeder", "facultative resprouter"), col=c("black", "red", "blue"), bty="n", pch=21)
plot(seq(0, 1, 0.1), k.prop.occ[,1], ylim=c(0,1), xlab="Conditional probability of a fire", ylab="Mean patch occupancy", bty="l")
points(seq(0, 1, 0.1), k.prop.occ[,2], col="red")
points(seq(0, 1, 0.1), k.prop.occ[,3], col="blue")
